#setwd()

library(foreign)
library(splines)
library(survival)

############################################################################
#Read Stata Data File                                                      #    
############################################################################
obscenity <- read.dta("condfrail_obscenity.dta")
summary(obscenity)

############################################################################
#Examine the Data                                                          #    
############################################################################

#Examine Data for Colorado
obscenity[41:48, 1:12]

#Examine Data for Florida
obscenity[65:72, 1:12]

#Number of Events in Data
eventcount <- table(obscenity$event)
eventcount

#Average Number of Events for Each State
avg.eventcount <- eventcount/max(obscenity$statenum)
avg.eventcount

#Frequency of Events in Each Strata
num.strata <- table(obscenity$evnum)
num.strata

#Maximum Number of Events for Each State
maxnum.state <- cbind(obscenity$statnam,obscenity$total.events)
maxnum.state

############################################################################
#Estimate a Conditional Cox Model in Gap Time                              #
############################################################################
cond.gap <- coxph(Surv(istart,ifinish,event)~fundam+south+collc+
                          legislativeideology+citideo+womleg+murderrate+
                          cluster(statenum)+strata(evnum),
                          method=c("efron"),data=obscenity)

###################################
#Obtain the Results from the Model#
###################################
summary(cond.gap)

#Collect the Estimated Beta's
cond.gap.beta <- matrix(cond.gap$coef,ncol=1)

#Collect the Robust S.E's
cond.gap.rse <- matrix(sqrt(diag(cond.gap$var)),ncol=1)

#################################
#Obtain Log-Likelihood for Model#
#################################
cond.gap.loglike2 <- cond.gap$loglik[[2]]
cond.gap.loglike2

#########################################
#Plot Cumulative Hazards for Each Strata#
#########################################
cond.gap.surv <- survfit(cond.gap) 
summary(cond.gap.surv)

plot(cond.gap.surv,fun="cumhaz",lty=c(1,2,3,4,5),col=c("black","black","black","black","black"),lwd=c(1,1,2,2,2),ylab="Cumulative Hazard",xlab="Time (Years)",ylim=c(0,4))
leg.txt <- c("Event 1","Event 2","Event 3","Event 4","Event 5")
legend(5,4.0,leg.txt,lty=c(1,2,3,4,5),col=c("black","black","black","black","black"),lwd=c(1,1,2,2,2))
temp <- c(1,2,3,4,5,6,7,8)
axis(1,temp)
box()

############################################################################
#Estimate a Cox Model in Elapsed Time with Shared Frailty from Gamma       #
#Distribution using the EM Algorithm                                       #
############################################################################
frailty.gamma.em <- coxph(Surv(estart,efinish,event)~fundam+south+collc+
                          legislativeideology+citideo+womleg+murderrate+ 
                          frailty(statename,distribution="gamma",method="em"),
			  method=c("efron"),data=obscenity)

#Estimate the same model without the frailty term to conduct an LR test for theta
nofrailty <- coxph(Surv(estart,efinish,event)~fundam+south+collc+
                          legislativeideology+citideo+womleg+murderrate,
                          method=c("efron"),data=obscenity)

###################################
#Obtain the Results from the Model#
###################################
summary(frailty.gamma.em)

#Collect the Estimated Beta's
frailty.gamma.em.beta <- matrix(frailty.gamma.em$coef,ncol=1)

#Collect the S.E's
frailty.gamma.em.se <- matrix(sqrt(diag(frailty.gamma.em$var)),ncol=1)

#################################
#Obtain Log-Likelihood for Model#
#################################
frailty.gamma.em.loglike2 <- frailty.gamma.em$loglik[[2]]
frailty.gamma.em.loglike2

################################################
#Obtain Estimated Variance of the Random Effect#
################################################
frailty.gamma.em.theta <- frailty.gamma.em$history[[1]]$theta
frailty.gamma.em.theta

###################
#LR test for theta#
###################
#Obtain I-likelihood
frailty.gamma.em.ill <- frailty.gamma.em$history[[1]]$c.loglik
frailty.gamma.em.ill

#Obtain likelihood from the no frailty model
nofrailty.loglike2 <- nofrailty$loglik[[2]] #Base likelihood for LR test
nofrailty.loglike2

#Calculate LR test for theta, report this test for whether theta is significant
frailty.gamma.em.lr <- 2*(341.58-340.95)
frailty.gamma.em.lr

#Calculate p-value, the test has 1 df 
frailty.gamma.em.lrtest <- 1-pchisq(frailty.gamma.em.lr,df=1)
frailty.gamma.em.lrtest

#################################################
#Calculate and Plot Profile Likelihood for theta#
#################################################
pllik.frailty.gamma.em <- double(100)
for(i in 1:100){
  a <- .01*i
  temp2 <- coxph(Surv(estart,efinish,event)~fundam+south+collc+
                          legislativeideology+citideo+womleg+murderrate+ 
                          frailty(statename,theta=a),
                          method=c("efron"),data=obscenity)
  pllik.frailty.gamma.em[i] <- temp2$history[[1]]$c.loglik
}

a <- seq(0.01,1,length=100)
temp <- frailty.gamma.em$history[[1]]$c.loglik - 3.84/2

#Calculate confidence interval using temp
frailty.gamma.em.ci <- cbind(a,pllik.frailty.gamma.em)

#Plot Profile Likelihood
plot(a, pllik.frailty.gamma.em, type='l', ylab="Log (Partial-Likelihood)",ylim=c(-352,-340),xlab=expression(paste(theta)))
abline(h=temp, lty=2)
points(frailty.gamma.em.theta,frailty.gamma.em.ill)
text(x=.23,y=-340.5,adj=c(.25,.25),c(expression(paste(hat(theta)))))
arrows(.22,-340.6,.185,-340.8,length=0.025)

######################################
#Obtain State Level Frailty Estimates#
######################################
frailties.gamma.em <- frailty.gamma.em$frail
summary(frailties.gamma.em)

#Obtain State Labels
state.abb <- c("AL","AK","AZ","AR","CA","CO","CT","DE","FL","GA","HI","ID","IL","IN","IA","KS","KY","LA","ME",
               "MD","MA","MI","MN","MS","MO","MT","NE","NV","NH","NJ","NM","NY","NC","ND","OH","OK","OR","PA",
               "RI","SC","SD","TN","TX","UT","VT","VA","WA","WV","WI","WY")

frailty.gamma.em.estimates.state <- cbind(state.abb,frailties.gamma.em)

######################################
#Plot Frailty Estimate for Each State#
######################################
plot(frailties.gamma.em,type='n',ylab="Frailty Estimates",xlab="State Number",ylim=c(-.4,.4))
text(frailties.gamma.em,state.abb)
abline(h=0,col="black",lty=2,lwd=1)
box()

########################
#Plot Cumulative Hazard#
########################
frailty.gamma.em.surv <- survfit(frailty.gamma.em,conf.type="none")
summary(frailty.gamma.em.surv)

plot(frailty.gamma.em.surv,fun="cumhaz",col="black",lwd=2, ylab="Cumulative Hazard",xlab="Time (Years)")
temp <- c(1,2,3,4,5,6,7,8)
axis(1,temp)
box()

############################################################################
#Estimate a Conditional Frailty Model in Gap Time with Shared Frailty from #
#Gamma Distribution using the EM Algorithm                                 #
############################################################################
condfrailty.gamma.em <- coxph(Surv(istart,ifinish,event)~fundam+south+collc+
                          legislativeideology+citideo+womleg+murderrate+
                          strata(evnum)+frailty(statename,distribution="gamma",method="em"),
                          method=c("efron"),data=obscenity)

#Estimate the same model without the frailty term to conduct an LR test for theta
condnofrailty <- coxph(Surv(istart,ifinish,event)~fundam+south+collc+
                          legislativeideology+citideo+womleg+murderrate+ 
                          strata(evnum),method=c("efron"),data=obscenity)

###################################
#Obtain the Results from the Model#
###################################
summary(condfrailty.gamma.em)

#Collect the Estimated Beta's
condfrailty.gamma.em.beta <- matrix(condfrailty.gamma.em$coef,ncol=1)

#Collect the S.E's
condfrailty.gamma.em.se <- matrix(sqrt(diag(condfrailty.gamma.em$var)),ncol=1)

#################################
#Obtain Log-Likelihood for Model#
#################################
condfrailty.gamma.em.loglike2 <- condfrailty.gamma.em$loglik[[2]]
condfrailty.gamma.em.loglike2

################################################
#Obtain Estimated Variance of the Random Effect#
################################################
condfrailty.gamma.em.theta <- condfrailty.gamma.em$history[[1]]$theta
condfrailty.gamma.em.theta

###################
#LR test for theta#
###################
#Obtain I-likelihood
condfrailty.gamma.em.ill <- condfrailty.gamma.em$history[[1]]$c.loglik #Report this one!
condfrailty.gamma.em.ill

#Obtain likelihood for the no frailty model
condnofrailty.loglike2 <- condnofrailty$loglik[[2]] #Base likelihood for LLR test
condnofrailty.loglike2

#Calculate LR test for theta, report this test for whether theta is significant
condfrailty.gamma.em.lr <- 2*(261.40-261.40)
condfrailty.gamma.em.lr

#Calculate p-value, it has 1 df
condfrailty.gamma.em.lrtest <- 1-pchisq(condfrailty.gamma.em.lr,df=1)
condfrailty.gamma.em.lrtest

#################################################
#Calculate and Plot Profile Likelihood for theta#
#################################################
pllik.condfrailty.gamma.em <- double(100)
for(i in 1:100){
  a <- 0.01*i
  temp2 <- coxph(Surv(istart,ifinish,event)~fundam+south+collc+
                          legislativeideology+citideo+womleg+murderrate+ 
                          strata(evnum)+frailty(statename,theta=a),method=c("efron"),data=obscenity)
  pllik.condfrailty.gamma.em[i] <- temp2$history[[1]]$c.loglik
}

a <- seq(.01,1,length=100)
temp <- condfrailty.gamma.em$history[[1]]$c.loglik - 3.84/2

#Calculate confidence interval using temp
condfrailty.gamma.em.ci <- cbind(a,pllik.condfrailty.gamma.em)

#Plot Profile Likelihood
plot(a, pllik.condfrailty.gamma.em, type='l', ylab="Log (Partial-Likelihood)",ylim=c(-272,-260),xlab=expression(paste(theta)))
abline(h=temp, lty=2)
points(condfrailty.gamma.em.theta, condfrailty.gamma.em.ill)
text(x=.05,y=-261,adj=c(.25,.25),c(expression(paste(hat(theta)))))
arrows(.04,-261,.01,-261.25,length=0.025)

######################################
#Obtain State Level Frailty Estimates#
######################################
condfrailties.gamma.em <- condfrailty.gamma.em$frail
summary(condfrailties.gamma.em)

condfrailty.gamma.em.estimates.state <- cbind(state.abb, condfrailties.gamma.em)

######################################
#Plot Frailty Estimate for Each State#
######################################
plot(condfrailties.gamma.em,type='n',ylab="Frailty Estimates",xlab="State Number",ylim=c(-.4,.4))
text(condfrailties.gamma.em,state.abb)
abline(h=0,col="black",lty=2,lwd=1)
box()

#########################################
#Plot Cumulative Hazards for Each Strata#
#########################################
condfrailty.gamma.em.surv <- survfit(condfrailty.gamma.em)
summary(condfrailty.gamma.em.surv)

plot(condfrailty.gamma.em.surv,fun="cumhaz",lty=c(1,2,3,4,5),col=c("black","black","black","black","black"),
     lwd=c(1,1,2,2,2),ylab="Cumulative Hazard",xlab="Time (Years)",ylim=c(0,4))
leg.txt <- c("Event 1","Event 2","Event 3","Event 4","Event 5")
legend(5,4,leg.txt,lty=c(1,2,3,4,5),col=c("black","black","black","black","black"),lwd=c(1,1,2,2,2))
temp <- c(1,2,3,4,5,6,7,8)
axis(1,temp)
box()

#######################################################################################
#Collect the Results the Conditional Gap Time, Frailty, and Conditional Frailty Models#                                                       
#######################################################################################

results.beta <- cbind(cond.gap.beta,frailty.gamma.em.beta,condfrailty.gamma.em.beta)
print.default(results.beta,2)

results.se <- cbind(cond.gap.rse,frailty.gamma.em.se,condfrailty.gamma.em.se)
print.default(results.se,2)

results.loglike <- cbind(cond.gap.loglike2,frailty.gamma.em.loglike2,condfrailty.gamma.em.loglike2)
print.default(results.loglike,6)

results.Ilike <- cbind(frailty.gamma.em.ill,condfrailty.gamma.em.ill)
print.default(results.Ilike,6)

results.theta <- cbind(frailty.gamma.em.theta,condfrailty.gamma.em.theta)
print.default(results.theta,2)

results.theta.lr <- cbind(frailty.gamma.em.lr,condfrailty.gamma.em.lr)
print.default(results.theta.lr,2)

results.theta.lrtest <- cbind(frailty.gamma.em.lrtest,condfrailty.gamma.em.lrtest)
print.default(results.theta.lrtest,2)
